Example SED analysis using PyXEM on DLS Jupyterhub

author: Mohsen Danaie

This notebook is intended to be run on Diamond's Jupyter hub. See instructions to connect here

This notebook covers the following tasks:

(a) Calibrations:

- setting up and loading files
- realspace pixel size calibration
- reciprocal space pixel size calibration
- correction of diffrcation roundness
- image / diffraction plane rotation calibration

(b) Example PyXEM functionalities
- Virtual dark-field imaging

(a) Calibrations

setting up the libraries needed

In [1]:
%matplotlib notebook
import numpy as np
import hyperspy.api as hs
import pyxem as pxm
import matplotlib.pyplot as plt
import os
import sys

from pyxem.libraries.calibration_library import CalibrationDataLibrary
from pyxem.generators.calibration_generator import CalibrationGenerator

import numpy as np
import matplotlib.pyplot as plt

from skimage import data
from skimage.registration import phase_cross_correlation
from skimage.registration._phase_cross_correlation import _upsampled_dft
from scipy.ndimage import fourier_shift

import math

import warnings
warnings.simplefilter(action='ignore')
plt.rcParams.update({'figure.max_open_warning': 0})
warnings.filterwarnings("ignore")

!git clone https://github.com/ePSIC-DLS/epsic_tools
current_path = os.getcwd()
sys.path.append(os.path.join(current_path, 'epsic_tools'))
import epsic_tools.api as epsic
/home/eha56862/.local/lib/python3.7/site-packages/pyUSID/viz/__init__.py:18: FutureWarning: Please use sidpy.viz.plot_utils instead of pyUSID.viz.plot_utils. pyUSID.plot_utils will be removed in a future release of pyUSID
  FutureWarning)
WARNING:silx.opencl.common:The module pyOpenCL has been imported but can't be used here
fatal: destination path 'epsic_tools' already exists and is not an empty directory.

loading files:

In [2]:
base_folder = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/' 
cals_path = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals'
mask_path = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/mask.hspy'

hdf_files = []
hspy_files = []
tiff_files = []
for dirname, dirnames, filenames in os.walk(cals_path):
    for filename in filenames:
        if os.path.splitext(filename)[1] == '.hdf5':
            hdf_files.append(os.path.join(dirname, filename))
        if os.path.splitext(filename)[1] == '.hspy':
            hspy_files.append(os.path.join(dirname, filename))
        if os.path.splitext(filename)[1] == '.tiff':
            tiff_files.append(os.path.join(dirname, filename))
for i, file in enumerate(hspy_files):
    print('[', i, ']', file.split('/')[-1])
[ 0 ] calibrated_dp.hspy
[ 1 ] mask.hspy
[ 2 ] vac_pattern.hspy
[ 3 ] au_xgrating_200kX_40cmCL_10umAp_subset_dp.hspy
[ 4 ] au_xgrating_200kX_40cmCL_10umAp_ibf.hspy
[ 5 ] defocus_scan_au_xgrating_200kX_40cmCL_10umAp_subset_dp.hspy
[ 6 ] defocus_scan_au_xgrating_200kX_40cmCL_10umAp_ibf.hspy
[ 7 ] defocus_scan_au_xgrating_200kX_40cmCL_10umAp_scan_array_125by127_diff_plane_515by515_flip.hspy
[ 8 ] defocus_scan_au_xgrating_200kX_40cmCL_10umAp_scan_array_125by127_diff_plane_515by515_flip_2.hspy
[ 9 ] binned_nav_au_xgrating_200kX_20cmCL_10umAp_scan_array_63by63_diff_plane_515by515__sumdiff.hspy
[ 10 ] au_xgrating_200kX_20cmCL_10umAp_ibf.hspy
[ 11 ] au_xgrating_200kX_20cmCL_10umAp_subset_dp.hspy
[ 12 ] binned_nav_au_xgrating_100kX_20cmCL_10umAp_scan_array_63by63_diff_plane_515by515__sumdiff.hspy
[ 13 ] au_xgrating_100kX_20cmCL_10umAp_ibf.hspy
[ 14 ] au_xgrating_100kX_20cmCL_10umAp_subset_dp.hspy
[ 15 ] binned_nav_au_xgrating_150kX_20cmCL_10umAp_scan_array_63by63_diff_plane_515by515__sumdiff.hspy
[ 16 ] au_xgrating_150kX_20cmCL_10umAp_subset_dp.hspy
[ 17 ] au_xgrating_150kX_20cmCL_10umAp_ibf.hspy
[ 18 ] au_xgrating_150kX_40cmCL_10umAp_subset_dp.hspy
[ 19 ] au_xgrating_150kX_40cmCL_10umAp_ibf.hspy
[ 20 ] au_xgrating_100kX_40cmCL_10umAp_ibf.hspy
[ 21 ] au_xgrating_100kX_40cmCL_10umAp_subset_dp.hspy

Calibrating real space pixel size

In [3]:
au_refImages = []

for file in hspy_files:
    if 'au_xgrating_100kX_20cmCL_10umAp_ibf' in file:
        au_refImages.append(file)
    if 'au_xgrating_150kX_20cmCL_10umAp_ibf' in file:
        au_refImages.append(file)
    if 'au_xgrating_200kX_20cmCL_10umAp_ibf' in file:
        au_refImages.append(file)
    

auRefIm = hs.load(au_refImages)

plt.figure()
plt.subplot(131, title = '100kX_20cmCL')
plt.imshow(auRefIm[1].data)
plt.axis('off')
plt.subplot(132, title = '150kX_20cmCL')
plt.imshow(auRefIm[2].data)
plt.axis('off')
plt.subplot(133, title = '200kX_20cmCL')
plt.imshow(auRefIm[0].data)
plt.axis('off')
Out[3]:
(-0.5, 254.5, 254.5, -0.5)

The pitch of this standard cross-grating sample is known (500 nm). Here we can measure and calibrate pixel sizes - interactively:

In [4]:
im= auRefIm[1]
im.plot()
line = hs.roi.Line2DROI(x1=68.7, y1=213.96, x2=189.9, y2=64.167, linewidth=10)
line.add_widget(im)      
trace = line.interactive(im)
trace.plot()
using interactive roi
In [5]:
MAG = [
    '100kX',
    '150kX',
    '200kX'
]
# pitch of the au cross-grating is 500 nm
nm_per_pixel = [
        np.mean([500 / (101-36.5),500 / (164.7-101)]),
        np.mean([500 / (109.46 - 14.91),500 / (205.97 - 109.46)]),
        np.mean([500 / (134.57 - 10.31),500 / (263.27 - 134.57)])
 ]
FOV = np.asarray(nm_per_pixel) * 256
#dictionary of MAG cal - MAG keys and FOV in nm values
MAG_cal_FOV = dict(zip(MAG, FOV))
print(MAG_cal_FOV)
# save to a file
np.save('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/MAG_cal_FOV_nm.npy', MAG_cal_FOV)
{'100kX': 1996.9576381536517, '150kX': 1340.034249786069, '200kX': 1012.3295878969468}

Calibrating reciprocal space pixel size

Reference diffraction data, only using 20 cm CL cases:

In [6]:
au_refDiff = []
for file in hspy_files:
    if '_sumdiff' in file:
        au_refDiff.append(file)
auRefDiff = hs.load(au_refDiff)


v_max = 15
plt.figure()
plt.subplot(131, title = '100kX_20cmCL')
plt.imshow(np.log(auRefDiff[2].data), vmax = v_max)
plt.axis('off')
plt.subplot(132, title = '150kX_20cmCL')
plt.imshow(np.log(auRefDiff[1].data), vmax = v_max)
plt.axis('off')
plt.subplot(133, title = '200kX_20cmCL')
plt.imshow(np.log(auRefDiff[2].data), vmax = v_max)
plt.axis('off')
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.
Out[6]:
(-0.5, 514.5, 514.5, -0.5)

It should not matter which one we choose to analyse here. So continuing with the first pattern:

In [7]:
dp = auRefDiff[0]

dp = pxm.signals.ElectronDiffraction2D(dp)
dp.plot(vmax = 8000)

Lets do a rough measurement of pixel size in reciprocal space - interactively:

In [8]:
dp.plot(vmax=8000)
line = hs.roi.Line2DROI(x1=186.251, y1=342.793, x2=307.451, y2=193.00, linewidth=5)
line.add_widget(dp)      

trace = line.interactive(dp)
trace.plot()
using interactive roi
In [9]:
# centre : 240, 275.5
# 111 ring diam in pixels: 92 - radius 46
# Au 111 reflection (A^-1): 1 / 2.355
diff_cal = (1/2.355) / 46
print(np.round(diff_cal, decimals=5), 'A^-1 to pix')
0.00923 A^-1 to pix

The above value is a rough initial estimate. We will refine this value later in the notebook.

Diffraction roundness correction

First we create a smaller dataset by skipping a number of probe positions, e.g. 10 here:

In [10]:
file = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_.hdf5'
dp_full = hs.load(file, lazy = True)
dp_skip10 = dp_full.inav[::10,::10]
# dp_skip10.save('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_skip10.hdf5')# dp_skip10= pxm.load_hspy('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_skip10.hdf5')
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.
Overwrite '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_skip10.hdf5' (y/n)?
n

Loading the mask and applying it to the skip10 data:

In [11]:
mask = hs.load(mask_path)
mask.plot()
dp_skip10_masked = dp_skip10 * mask
dp_skip10_masked = pxm.signals.ElectronDiffraction2D(dp_skip10_masked)
dp_skip10_masked_sum = dp_skip10_masked.sum()
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.

We need to have the central beam at the centre of the diffraction patterns to continue with this analysis. Here we find the centre of the pattern on the sum diffraction pattern and crop around this position:

In [12]:
centre = pxm.utils.expt_utils.find_beam_center_interpolate(dp_skip10_masked_sum, sigma=5, upsample_factor=4, kind='linear')
print('Coordinates of the centre: ', centre)
det_size = dp_skip10_masked_sum.axes_manager[0].size
crop_size = min(int(det_size - centre[0]),int(det_size - centre[1]), centre[0], centre[1])
print('croped size of the patterns: ', int(2 * crop_size))
dp_skip10_masked_crop = dp_skip10_masked.isig[int(centre[0]-crop_size):int(centre[0]+crop_size), int(centre[1]-crop_size):int(centre[1]+crop_size)]
Coordinates of the centre:  [276.25 241.25]
croped size of the patterns:  476

Now we can correct for shifts in the diffraction pattern by centring the direct beam:

In [13]:
dp_skip10_masked_crop.center_direct_beam(method='interpolate', sigma=5, upsample_factor=4, kind='linear')
dp_skip10_masked_crop.axes_manager[2].offset = 0
dp_skip10_masked_crop.axes_manager[3].offset = 0
[########################################] | 100% Completed |  1.7s
In [14]:
dp_skip10_masked_crop.plot()

saving a centred version:

After the above correction, we can sum all the diffraction patterns and apply the reciprocal space calibration:

In [15]:
dp_sum = dp_skip10_masked_crop.sum()
dp_sum = pxm.signals.ElectronDiffraction2D(dp_sum)
dp_sum.set_experimental_parameters(beam_energy=300, camera_length=20.)
dp_sum.set_diffraction_calibration(calibration=diff_cal)
cal_lib = CalibrationDataLibrary(au_x_grating_dp=dp_sum)
# cal = CalibrationGenerator(calibration_data=cal_lib)
cal = CalibrationGenerator(diffraction_pattern=dp_sum)
dp_sum.plot(vmax = 300, cmap='magma_r')

Using this pattern we can correct for diffraction roundness:

In [16]:
cal.get_elliptical_distortion(mask_radius=10,
                              scale=diff_cal, amplitude=80, direct_beam_amplitude=4e3,
                              asymmetry=0.9,spread=2)
corr_matrix = cal.get_correction_matrix()
print(corr_matrix)
residuals = cal.get_distortion_residuals(mask_radius=10, spread=2)
residuals.plot(cmap='RdBu', vmax=0.005)
[[0.9920649  0.01695637 0.        ]
 [0.01695637 0.96376626 0.        ]
 [0.         0.         1.        ]]
In [17]:
cal.plot_corrected_diffraction_pattern(vmax= 80)

Now we can refine the reciprocal space pixel size

In [18]:
centre = (int((dp_sum.data.shape[0]/2)), int((dp_sum.data.shape[1]/2)))
profile = epsic.radial_profile.radial_profile(dp_sum.data, centre)
fig, ax = plt.subplots(1,1)
ax.plot(profile)
ax.set_xlim([0, 50])
ax.set_ylim([0, 200])
ax.annotate('au_111', xy=(44, 95),  xycoords='data',
            xytext=(0.8, 0.95), textcoords='axes fraction',
            arrowprops=dict(facecolor='black', shrink=0.05),
            horizontalalignment='right', verticalalignment='top',
            )
Out[18]:
Text(0.8, 0.95, 'au_111')
In [19]:
diff_cal_refined = (1/2.355) / 44.4
print('Initial estimate:', np.round(diff_cal,5), 'A^-1 to pix')
print('Refined value: ', np.round(diff_cal_refined,5), 'A^-1 to pix')
Initial estimate: 0.00923 A^-1 to pix
Refined value:  0.00956 A^-1 to pix

Scan rotation calibration

This calibration aims to correct for the rotation angle between the scan coil fast and slow directions and the X/Y axis of the detector. The data we load here is a 4DSTEM dataset of a largely defocused probe (i.e. Ronchigram). We then take two lines corresponding to the slow and fast scan directions.

In [20]:
scan_data = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 115307/from_mg25142-9/20201214_153300.hdf5'
scan = hs.load(scan_data, lazy = True)
fast_dir = scan.inav[1:,60]
slow_dir = scan.inav[60,1:]
fast_dir.compute()
slow_dir.compute()
[########################################] | 100% Completed |  1.3s
[########################################] | 100% Completed |  1.8s
In [21]:
slow_dir.plot()

Cropping to the area within the bright-field disc, removing the cross, and applying a Gaussian blur to better track the motion:

In [22]:
index = [255,256,257,258,259]
fast_dir_no_cross = np.delete(fast_dir.data, index, axis = 1)
fast_dir_no_cross = np.delete(fast_dir_no_cross, index, axis = 2)
slow_dir_no_cross = np.delete(slow_dir.data, index, axis = 1)
slow_dir_no_cross = np.delete(slow_dir_no_cross, index, axis = 2)
fast_dir_no_cross = hs.signals.Signal2D(fast_dir_no_cross)
slow_dir_no_cross = hs.signals.Signal2D(slow_dir_no_cross)
fast_dir_crop = fast_dir_no_cross.isig[132:429, 108:415]
slow_dir_crop = slow_dir_no_cross.isig[132:429, 108:415]
from scipy.ndimage import gaussian_filter
fast_dir_crop_gb = gaussian_filter(fast_dir_crop.data, sigma = 2)
fast_dir_crop_gb = hs.signals.Signal2D(fast_dir_crop_gb)
slow_dir_crop_gb = gaussian_filter(slow_dir_crop.data, sigma = 2)
slow_dir_crop_gb = hs.signals.Signal2D(slow_dir_crop_gb)
slow_dir_crop_gb.plot()

estimate the shifts in the image:

In [23]:
shifts_slow = slow_dir_crop_gb.inav[:20].estimate_shift2D('current', roi = (132,292,150,300), normalize_corr=True)
shifts_fast = fast_dir_crop_gb.inav[:20].estimate_shift2D('current', roi = (132,292,150,300), normalize_corr=True)
In [24]:
shifts_slow = np.asarray(shifts_slow)
print(np.sum(shifts_slow, axis=0))
shifts_slow_sum = np.sum(shifts_slow, axis=0)
shifts_fast = np.asarray(shifts_fast)
print(np.sum(shifts_fast, axis=0))
shifts_fast_sum = np.sum(shifts_fast, axis=0)
[ 42 160]
[-170   30]
In [25]:
plt.figure()
plt.plot([0, shifts_slow_sum[0]], [0, shifts_slow_sum[1]], '-r', label = 'slow_dir')
plt.plot([0, shifts_fast_sum[0]], [0, shifts_fast_sum[1]], '-b', label = 'fast_dir') 
plt.axis('square')
plt.legend()
Out[25]:
<matplotlib.legend.Legend at 0x7f943935d3d0>
In [26]:
import scipy
rot_angle = (abs(math.degrees(math.atan(shifts_fast_sum[1]/shifts_fast_sum[0])))
            + 180 - abs(math.degrees(math.atan(shifts_fast_sum[1]/shifts_fast_sum[0])))
            - abs(math.degrees(math.atan(shifts_slow_sum[1]/shifts_slow_sum[0]))))
print(rot_angle)
slow_dir_crop_gb_rot = scipy.ndimage.rotate(slow_dir_crop_gb.data, angle=rot_angle, axes = (1,2))
slow_dir_crop_gb_rot = hs.signals.Signal2D(slow_dir_crop_gb_rot)
slow_dir_crop_gb_rot.plot()
104.70830389968275
In [27]:
fast_dir_crop_gb_rot = scipy.ndimage.rotate(fast_dir_crop_gb.data, angle=rot_angle, axes = (1,2))
fast_dir_crop_gb_rot = hs.signals.Signal2D(fast_dir_crop_gb_rot)
fast_dir_crop_gb_rot.plot()

(b) Example PyXem functionalities

Forming virtual dark-field images

In [28]:
test_data = hs.load(hdf_files[0])
test_data = pxm.signals.ElectronDiffraction2D(test_data)
test_data.plot(vmin=0.,vmax=10)
WARNING:hyperspy.io:`signal_type='STEM'` not understood. See `hs.print_known_signal_types()` for a list of known signal types, and the developer guide for details on how to add new signal_types.

The green circular region of interest above (this appears after running the cell below), defines a virtual aperture. The sum intensity of this region then is displayed as a virtual dark-field image below. By changing the size and location of the aperture size, the VDF image is updated interactively.

In [29]:
roi = hs.roi.CircleROI(cx=313, cy=308, r=10, r_inner=0)
test_data.plot_integrated_intensity(roi)
using interactive roi